#include<iostream>
using namespace std;
#include<iomanip>
#include<cmath>

#include"library.h"

/*
 * SelfAdapt_trapzoid uses self adapted algorithm to calculate integration of func() over [a,b].
 * It automatically breaks [a,b] into unequal pieces, according to the local situaion of func().
 * the error is controled to be less than epsilon
 */
double SelfAdapt_trapzoid(double a, double b, double (*func)(double), double epsilon){

        double y1 = integ_trapzoid(a, (a+b)/2, 1, func);
        double y2 = integ_trapzoid((a+b)/2, b, 1, func);
        double y3 = integ_trapzoid(a, b, 1, func);
        if( fabs(y1+y2-y3) < 3*epsilon ){// 3, instead of 15 here, cause of similar derivations of SelfAdapt_Simpson
                return y3;
        }
        else{
                return SelfAdapt_trapzoid(a, (a+b)/2, func, epsilon/2)
                        +
                        SelfAdapt_trapzoid((a+b)/2, b, func, epsilon/2);
        }
}

#define exact 1.15288726707017

double y(double x){
        if(fabs(x)<1E-9){
                cout<<"invalid input for y(x): |x|<1E-9"<<endl;
                exit(1);
        }
        return sin(1/x/x);
}

int main(){

        double epsilon = 1E-10;
        double a = 0.1, b = 10;
        double z1 = SelfAdapt_Simpson(a,b,y,epsilon);
        cout<<setprecision(15)<<"SelfAdapt_Simpson: error="<<z1-exact<<"\t";

        double z2;
        cout<<"Simpson:"<<endl;
        for(int i=0;i<8;i++){
                z2 = Simpson(a,b,pow(10,i),y);
                cout<<"h="<<1/pow(10,i)<<"\t\t";
                cout<<pow(10,i)<<" grids\t";
                cout<<"error="<<z2-exact<<endl;
        }

        return 0;
}
